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Abstract 

Background: Statistical process control (SPC), an industrial sphere initiative, has recently been applied in health 
care and public health surveillance. SPC methods assume independent observations and process autocorrelation 
has been associated with increase in false alarm frequency. 

Methods: Monthly mean raw mortality (at hospital discharge) time series, 1995-2009, at the individual Intensive Care 
unit (ICU) level, were generated from the Australia and New Zealand Intensive Care Society adult patient database. 
Evidence for series (i) autocorrelation and seasonality was demonstrated using (partial)-autocorrelation ((P)ACF) function 
displays and classical series decomposition and (ii) "in-control" status was sought using risk-adjusted (RA) exponentially 
weighted moving average (EWMA) control limits (3 sigma). Risk adjustment was achieved using a random coefficient 
(intercept as ICU site and slope as APACHE III score) logistic regression model, generating an expected mortality series. 
Application of time-series to an exemplar complete ICU series (1995-(end)2009) was via Box-Jenkins methodology: 
autoregressive moving average (ARMA) and (G)ARCH ((Generalised) Autoregressive Conditional Heteroscedasticity) 
models, the latter addressing volatility of the series variance. 

Results: The overall data set, 1995-2009, consisted of 491324 records from 137 ICU sites; average raw mortality was 
14.07%; average(SD) raw and expected mortalities ranged from 0.012(0.1 13) and 0.013(0.045) to 0.296(0.457) and 0.278 
(0.247) respectively. For the raw mortality series: 71 sites had continuous data for assessment up to or beyond lag 40 and 
35% had autocorrelation through to lag 40 ; and of 36 sites with continuous data for > 72 months, all demonstrated 
marked seasonality. Similar numbers and percentages were seen with the expected series. Out-of-control signalling was 
evident for the raw mortality series with respect to RA-EWMA control limits; a seasonal ARMA model, with GARCH 
effects, displayed white-noise residuals which were in-control with respect to EWMA control limits and one-step 
prediction error limits (3SE).The expected series was modelled with a multiplicative seasonal autoregressive model. 

Conclusions: The data generating process of monthly raw mortality series at the ICU level displayed autocorrelation, 
seasonality and volatility. False-positive signalling of the raw mortality series was evident with respect to RA-EWMA 
control limits. A time series approach using residual control charts resolved these issues. 

Keywords: Statistical process control, Time series, Autocorrelation, Seasonality, Volatility, Exponentially weighted moving 
average smoothing, Autoregressive moving average models, GARCH models 
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Background 

Statistical process control (SPC), deriving from Shewart s 
work in 1920-30 and in the 1950s with Demings refine- 
ments [1], has been more recently applied in health care 
and public health surveillance [2], generating consider- 
able interest in the general [3-5] and specialist medical 
literature [6-10]; and has been subject to a detailed ex- 
position from a "quality-in-medicine" perspective [11]. 
Important statistical principles underlying SPC or 
control-chart methodology are those of the monitored 
process being "in control" and subject to the independ- 
ence of observations [12]. The presence and impact 
(possible increase in frequency of false alarms) of 
process autocorrelation in industrial/engineering series 
have long been documented [13-16]. Somewhat surpris- 
ingly, little formal attention has been directed to this 
problem in the bio-medical literature [17,18], one review 
suggesting that there was "...limited advice on how to 
manage [autocorrelation]..." [5]. 

We have previously drawn attention to the data- 
generating mechanisms of overall monthly mortality 
series, at the aggregate level, from a bi-national intensive- 
care (ICU) database, where persistent autocorrelation 
(to lag 2 4) was evident in a seasonal ARIMA (auto- 
regressive integrated moving average) model of the 
mortality series [19]. We now extend this study to 
further characterise the data generating process of 
mortality series at the individual ICU level and the 
impact of autocorrelation upon (i) mortality monitor- 
ing using EWMA (exponentially weighted moving 
average) control charts and (ii) time-series modelling 
of the data process using residual control charts. 

Methods 

As previously described [19,20], the ANZICS (Australian 
and New Zealand Intensive Care Society) adult patient 
database [21] was utilised to define an appropriate pa- 
tient set, 1995-(end)2009. Physiological variables col- 
lected in accordance with the requirements of the 
APACHE III algorithm [22,23] were the worst in the 
first 24 hours after ICU (intensive care unit) admis- 
sion, and all first ICU admissions to a particular hos- 
pital for the period 1995-2009 were selected. Records 
were used only when all three components of the 
Glasgow Coma Score (GCS) were provided; records 
for which all physiologic variables were missing were 
excluded, and for the remaining records, missing variables 
were replaced with the normal range and weighted ac- 
cordingly. The mortality endpoint was at hospital dis- 
charge. Exclusions: unknown hospital outcome; patients 
with an ICU length of stay < 4 hours, and patients 
aged < 16 years of age. Access to the data was 
granted by the ANZICS Database Management Com- 
mittee in accordance with standing protocols; local hospital 



(The Queen Elizabeth Hospital) Ethics of Research Com- 
mittee approval was waived. 

Statistical analysis 

(i) Monthly raw (risk-unadjusted) and risk- adjusted 
(RA) mortality time series at the individual ICU 
were generated. Risk adjustment was undertaken, 
generating the "expected" series, using a random 
coefficient logistic model (intercept as ICU site and 
"slope" as (centred) APACHE III score; 
unstructured covariance using adaptive quadrature, 
estimated via the Stata™ module "xtmelogit" [24]), 
as previously described in detail [20], and extended 
to both ventilated and non-ventilated patients. No 
formal adjustment for potential seasonality 
(trigonometric seasonality using sine/cosine 
functions or monthly dummy variables) was 
undertaken. Individual ICUs were allocated an 
identifier based upon a random number sequence. 

(ii) Graphical inspection of the mortality series and 
formal testing of normality to confirm that the 
"...distributions of... (observed) and... (predicted) 
[series] ... were sufficiently similar and are 
robustly normal and symmetrical..." [25]. 
Classical seasonal decomposition [26] was 
undertaken using the "decompose" module in R 
statistical software (Version 15.2 [27]). 
Autocorrelation plots (scatterplot grid of series 
versus lagged values) were performed via the R 
user-written module "lagl.plot" [28]. 

(iii) Generation of EWMA charts with confidence limits. 

a. assuming iid (independent and identically 
distributed) observations, the EWMA statistic fe) 
is defined as: Xx t + (1 - A)^_i and the variance 
(a 2 ) as a\ = a x 2 (£) [l-(l-\) 2i ] , where 0 < X < 1 
is a constant (smoothing parameter) [29] . 

b. For the variance of (non-stationary) auto-correlated 
series, we followed Montgomery & Mastrangelo 
[15]: division of the sum of squared (prediction) 
errors for optimal \ by n; leading to the plotting of 
a moving centre-line EWMA control chart [12]. 

c. Default values ("optimal") in Stata™ statistical 
software for \ were chosen to minimize the in- 
sample sum-of-squares forecast errors [30], a 
method also recommended by Montgomery and 
Mastrangelo [15]; albeit small values of X may 
inhibit the detection of large sudden process 
shifts; the "inertia" phenomenon [31]. 

d. Average run length (ARL): that is, the average 
number of "points", when the data-generating 
process is in fact in-control, plotted before out- 
of-control is declared (ARL 0 ). For instance, with 
iid observations and a Shewhart control-chart 
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with three sigma limits, ARL 0 = 1/^=1/0.0027=370 
(where p is the probability that any point 
exceeds the control limits [32,33], when the 
data-generating process is in fact in-control). 
Under the iid assumption, for various 
mortality series and values of \, scenario based 
increments of the (mean of the) underlying 
series were computed using Statgraphics® 
Centurion XVI statistical software [34]. 
e. Using conventional SPC methods, EWMA 
control limits (at 3 sigma) were applied to the 
raw mortality series using the expected series as 
reference process; that is, RA control limits were 
generated. 

(iv) Establishment of time-series models at the individual 
ICU level was based upon classic Box-Jenkins 
methodology (autoregressive moving average 
(ARMA) models) with investigation of (G)ARCH 
((Generalised) Autoregressive Conditional 
Heteroscedasticity) effects [35,36], as previously 
described [19]. 

a. A stationary time series {x t ; t = 0, ± 1, ± 2, ...} 
has an autoregressive moving average (ARMA 
(p,q)) structure: x t = (p 1 x t _ 1 + ... (ppX t _ p + % + 
0iG)f_i + ... 0 q o) t _ q where fa, fa, (j) p are the 
"autoregressive" (AR) coefficients relating the 
value of x at time t to its past p values, and 
0i, 6 2 , 6 q are the "moving average" (MA) 
coefficients, relating the current "white- 
noise",%, to its past q values and co t ~N(0, of,). 
If x t has a non-zero mean (pt), then a constant 
a = ^(1 - fa - ... - (p p ) is introduced into the 
structure. An integrated series accumulates 
(some) past effects and is therefore non-stationary. 
A series is integrated, say, of order 1 (1(1)) if the 
changes (or differences: Ax t = x t -x t _ 1 ) of the series 
generate stationarity (1(0)), leading to the 
expanded ARIMA model (ARlMA(p,d,q)), 
where d is the degree of differencing [37]. 
This being said, careful attention was directed 
to the question of trend versus difference 
stationarity [38], especially in medical series 
where, as opposed to stochastic random walks, 
"deterministic" trends may be present. [39,40]. 

b. Model diagnostics: the use of auto- (ACF) and 
partial-autocorrelation (PACF) function 
displays, testing for the presence of a unit-root 
(ADF (augmented Dickey-Fuller) and DF-GLS 
(modified Dickey-Fuller t test) tests [30] and 
variants), residual white-noise (Bartletts 
periodogram-based- and Portmanteau (Q)-test) 
and seasonality were undertaken after 
Shumway & Stoffer [41] and as previously 
described [19]. 



c. Volatility of the (squared) residuals (e) of the 
mean equation (conditional heteroscedasticity 
[42]) was checked using the PAC of the squared 
residuals and the user-written Stata™ "armadiag" 
module [43]; that is, ARCH and GARCH effects 
((Generalised) Autoregressive Conditional 
Heteroscedasticity of the error variance process). 
For an ARCH model, the mean equation is y t = 
Xtfi + e t and the variance equation of = y 0 + y 1 
e^_ 1 + y 2 £ t-2 + wnere £t~N(0, of), are the 
squared residuals (innovations) and ji are the 
ARCH parameters; the conditional variance is thus 
modelled as an AR process. A GARCH(m,/c) model 
includes lagged values of the conditional variance 

0? = y 0 + Yie?_l + Y2 £ t 2 -2 + - + Ym £ t-m + 
&l G t-l + &2 G t-2 + + Sk^t-k)' wnere $i are tne 

GARCH parameters (an ARMA process) [19,44]. 
Exploration of different error term distributions 
(normal, t and generalised error) was also 
undertaken [30]. 

d. Under the conditions of an appropriately 
specified time-series model, the behaviour of 
the residuals was investigated, after Alwan and 
Roberts [45], on the basis that a shift in the 
mean of a time series is transmitted to the 
residuals [46]. 

i. As residuals are assumed to be independent 
(white-noise: a sequence of iid random 
variables with finite mean and variance, all 
ACFs being [close to] zero [47]), standard 
control chart methods were used to generate 
residual-EWMA charts [33]. Thus, 
determination of the residual-EWMA 
smoothing parameter (X) was based upon 
methods for independent observations. 

ii. Control limits were also determined using 
standard errors (3x) of the one-step-ahead 
forecasts [45]. 

e. Model selection was guided by penalized 
information criteria (Akaike (AIC) and Bayesian 
(BIC) information criteria) [48]. 

f. Formal exegesis proceeded using a single 
exemplar complete ICU series (1995-(end)2009). 

(v) Graphical displays: line-graphs of series were 
produced for appropriate illustration of relevant 
stages of analysis 

a. Line graph(s) of the raw series were produced 
with 3*SE control limits of the expected series. 

b. EWMA control limits (including residual control 
charts) were generated using default values of 
"optimal exponential coefficient" in Stata™ 
statistical software [49]. 

c. Values of X for scenario based increments (say, 
5% or 10%) of target mean were calculated using 
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the SPC module of Statgraphics® statistical 
software [34] and appropriate 3*SE control limits 
of the expected series as in (a) above or EWMA 
line graphs were produced as in (b) above. 

Results 

The overall data set, 1995-2009, consisted of 491324 re- 
cords from 137 ICU sites; mean (hospital) mortality was 
14.07%. The random coefficient logistic regression model 
(Hosmer-Lemehsow statistic 62.97, ROC area under the 
curve 0.89) generated an overall predicted mortality prob- 
ability of 0.1407 (SD 0.0202, range 0.00004-0.993). Over 
the 137 sites mean raw and expected (RA-) mortalities 
ranged from 0.012(0.113) and 0.013(0.045) to 0.296(0.457) 
and 0.278(0.247) respectively. 

Of the raw mortality series from the 137 ICUs, 71 had 
continuous monthly data (excluding missing values or zero 
monthly mortality) for assessment up to or beyond lag 40 . 
For 25 of these series (35%), there was a significant Q test 
(null hypothesis being that the series is white noise) and 
autocorrelation through to lag4 0 . Thirty six had continuous 
monthly data (excluding missing values) for > 72 months; 
all series demonstrated marked seasonality and 30 demon- 
strated an obvious trend decline in mortality. Of the 
expected mortality series, 72 had appropriately assessable 
data to lag4o and in 46 (64%) there was a significant Q test 
and autocorrelation through to lag 40 . Similarly, in the same 
36 series with continuous (raw) monthly data for > 
72 months, all expected mortality series demonstrated 



marked seasonality and 30 demonstrated an obvious trend 
decline in mortality. 

Data from site "4" over 1995-2009 was used to gen- 
erate an exemplar mortality time series. The mean 
raw mortality was 0.139(0.047) with skewness 0.216 
and kurtosis 2.53; and the expected mortality was 
0.138(0.028) with skewness 0.361 and kurtosis 3.47. 
The Shapiro-Wilk normality test was not rejected for 
either series (P =0.23 for both series). Kernel density 
estimates of raw and expected mortality are seen in 
Figure 1 (upper panel), with obvious difference in the 
degree of kurtosis between the two series. Time series 
plots, 1995-2009, for raw and expected (RA-)mortality 
are seen in the lower panel; a gradual time-decline in 
mortality for both series is evident. Additive seasonal 
decomposition of both series is seen in Figure 2, re- 
vealing marked seasonality and a trend decline in 
mortality. Autocorrelation plots are seen in Figure 3, 
showing correlation (positive and negative) decreasing 
variably with increase in lag in both series. 

Figure 4 displays a plot of raw mortality series with 
control limits as 3SE of expected mortality (upper 
panel) and a scenario based mortality increment of 5% 
(5% false positive rate and desired ARL= 6 months) with 
control limits as 3SE of expected mortality. Frequent 
signalling is seen in both panel-plots. Figure 5 shows a 
plot of the raw mortality series with (fixed) EWMA 3 SE 
control limits derived from a projected 5% (upper panel) 
and 10% (lower panel) increment in expected mortality, 



Monthly mortality: raw & expected 




.125 .15 .175 

Mortality 



Raw (hospital) mortality 
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Raw & expected monthly mortality series: 1995-2009 
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Raw (hospital) mortality 
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Figure 1 Upper panel: Kernel density estimates of raw and expected mortality. Lower panel: Time series (1995-2009) of mean monthly 
raw and expected mortality. 



Moran and Solomon BMC Medical Research Methodology 2013, 13:66 
http://www.biomedcentral.com/1471-2288/13/66 



Page 5 of 12 



Raw mortality 

Decomposition of additive time series 




Time 



Expected mortality 

Decomposition of additive time series 
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Figure 2 Decomposition plot of raw (upper panel) and 
expected (lower panel) mortality series. 



assuming: an in-control ARL of 370, mean (expected) 
mortality 0.1381(0.0276) and target mean (expected) 
mortality of 0.145 (5% increment) and 0.152 (10% incre- 
ment), for an EWMA \ = 0.02 and 0.05, respectively 
(calculations preformed in Stagraphics®). For both 5% 
and 10% projected increments of expected mortality, 
the raw mortality series signalled frequently, mainly in 
the early periods. Figure 6 shows the same scenarios 
with a time-varying variance EWMA control chart; 
again, there was frequent signalling of the raw mortality 
series. 



The autocorrelation evident in the raw and expected 
mortality series suggested a formal time series approach 
to SPC: 

(i) Raw mortality: both the DF-GLS and ADF tests 
(with trend) rejected the null-hypothesis of 
presence of a unit-root and the series was de- 
trended using linear regression (raw mortality 
against time) and the residuals (also not 
evidencing a unit-root) of the linear regression 
model were used for subsequent formal analysis. 
The de-trended series from the raw mortality 
displayed seasonality but, not surprisingly, no 
trend decline (graphics not shown). An initial 
additive seasonal ARMA model satisfied 
conventional diagnostic requirements, but 
displayed ARCH effects. Of the (G)ARCH models 
assessed, the most parsimonious was a simple 
[ARCH-lagi, GARCH-lagi] model (Table 1). 
Although the individual GARCH term was 
nominally non-significant, there was a highly 
significant (P=0.0001) test of joint significance 

of the ARCH and GARCH parameters. There 
was no advantage of either t or general error 
distribution in the development of the (G)ARCH 
models. 

(ii) Expected mortality: trend stationarity was 
demonstrated by rejection of the null-hypothesis of 
existence of a unit-root by the DF-GLS and ADF 
tests (with the trend option) and de-trending (linear 
regression of expected mortality against time) 
yielded residuals (also not evidencing a unit-root) 
for subsequent formal analysis. A simple 
(multiplicative) seasonal autoregressive model was 
generated with no evidence of ARCH effects 
(Table 1). Although an ARMA(1,1) model satisfied 
model diagnostic tests, the multiplicative seasonal 
AR model was favoured on clinical grounds. 

Both the GARCH and ARMA models were considered 
parsimonious and the de-trended signals for each model 
were within 3SE limits of respective model predictions 
(Figure 7). The residuals from both the formal GARCH 
and ARMA models (mean: 0(0.0423) and 0(0.0257) re- 
spectively) satisfied multiple criteria of Gaussian white- 
noise and were within residual-EWMA control limits 
(default values of "optimal exponential coefficient" in 
Stata™ statistical software; 3SE control limits; \ = 0.0001 
for both series; Figure 8). To address any potential iner- 
tial problems consequent upon the small A, control 
limits were also established for projected 1 (X=0.16), 2 
(\=0.42) and 3 (\=0.71) SD increments of the mean of 
the GARCH residuals; the latter were within these con- 
trol limits (Figure 9). 
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Figure 3 Lagplot of series versus lagged values (to lag 24 ); upper panel, raw mortality; lower panel, expected mortality. 
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Raw monthly mortality series; 3SE expected CI: 1995-2010 
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Figure 4 Upper panel: raw mortality series with 3SE control limits of expected mortality; lower panel: EWMA (X=0.51) of raw mortality series 
for anticipated 5% increase in raw mortality, 5% fale positive rate and desired ARL=6 months with 3SE control limits of expected mortality. 



Discussion due to Simpsons paradox. We thus concur with the 

The current analysis of monthly mortality series confirms findings of Alwan [13,50] and Bisgaard and Kulahci [51], 

the existence of autocorrelation and seasonality in both the who documented the pervasiveness of autocorrelation in a 

raw and expected series at the individual ICU level, variety of series, industrial and non-industrial. We also 

avoiding any potential confounding at the aggregate level established that out-of-control signalling of the raw 



EWMA A, = 0.02 for expected mortality (5% increase), fixed control limits for raw mortality 
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EWMA X = 0.05 for expected mortality (10% increase), fixed control limits for raw mortality 




Calendar time: monthly 

Figure 5 Raw mortality series with EWMA fixed 3SE control limits: upper control limit (red line), lower control limit (green line), signal 
(navy line); for anticipated 5% (upper panel) and 10% (lower panel) increments in expected mortality. 
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EWMA X = 0.02 for expected mortality (5% increase), time-varying control limits for raw mortality 
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EWMA X = 0.05 for expected mortality (10% increase), time-varying control limits for raw mortality 
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Figure 6 Raw mortality series with EWMA time-varying 3SE control limits: upper control limit (red line), lower control limit (green line), 
signal (navy line); for anticipated 5% (upper panel) and 10% (lower panel) increments in expected mortality. 



mortality series with respect to both 3 standard error risk- 
adjusted and RA-EWMA control limits was not evident 
with analysis of the residuals from the GARCH time series 
model Thus the identification of (G)ARCH processes is an 
important issue for SPC [35]. 

As our focus was directed to an understanding of the 
underlying data-generating process [45] and the perform- 
ance of the RA-EWMA control limits under conditions of 
autocorrelation, we deemed it appropriate to also subject 
the expected series, from which the control limits for the 
raw mortality series were established, to formal time-series 
estimation. Not surprisingly, as the underlying mortality 
estimates from a random coefficient model are obligatorily 
"smoothed" (see also Figure 1), no ARCH effects, 
representing "volatility", were demonstrated and a relatively 
simple seasonal autoregressive model was established 
(Table 1). As the EWMA is based upon an ARIMA(0,1,1), 
that is an integrated moving average process [45,52], it has 
been applied to autocorrelated data [15], although the ma- 
jority of studies have used relatively simple non-seasonal 
autoregressive models (AR(1) or AR(2)) with fixed \ 
(usually 0.2, which is the default for the SPC model of 
Statgraphics software). Residual-EWMA charts, in the 
context of time series modelling, would appear to be more 
robust than EWMA applied to the original (autocorrelated) 
data [53,54]. Reynolds and Lu have recommended that 
under autocorrelation "...traditional control chart method- 
ology should not be applied without modification..." [55] 



Table 1 Parameters for GARCH (estimated from raw 
mortality series; de-trended linear model residuals) and 
ARMA (estimated from the expected mortality series; 
de-trended linear model residuals) models 



Model 

Raw mortality 



Coefficient P 



UL_95%CI LL_95%CI 



ARMA 

Autoregressive parameters 

L24 0.174 0.016 0.033 0.316 

Moving average parameters 

L1 -0.148 0.021 -0.274 -0.022 

L15 0.265 0.000 0.128 0.402 

L17 -0.203 0.003 -0.340 -0.067 

ARCH 

L1 0.014 0.464 -0.024 0.052 

GARCH 

L1 0.996 0.000 0.931 1.061 

Expected mortality 

ARMA 

Autoregressive parameters 

L1 0.145 0.045 0.003 0.288 

ARMA12 

Autoregressive parameters 

L1 0.132 0.075 -0.031 0.277 

L lag. 
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Figure 7 De-trended series (navy line) generating the GARCH (upper panel) and ARMA (lower panel) models with one-step-ahead 
forecast control limits (3SE); upper control limit (red line), lower control limit (green line). 



and Human et al. have recently sounded a cautionary note 
about the robustness of the conventional EWMA [56]. 

For the classical SPC model, a process is in control if 
the mean and standard deviation estimate remain within 
prescribed control limits [57], usually three-sigma; that 



is, for a normally distributed series, 99.7% of observa- 
tions should lie within the limits [58] and there a prob- 
ability of 0.0027 that any point exceeds the control limit 
[32,59]. However this definition does not necessarily entail 
the formal time-series notion of stationarity (strict or weak), 



Detrended raw series vs GARCH predictions, 3SE CL 
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Figure 8 "Optimal" residual-EWMA control chart (3 SE control limits): upper control limit (red line), lower control limit (green line), 
residuals (navy line); for GARCH (upper panel) and ARMA (lower panel) model residuals, respectively. 
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where the requirement for stationarity is that the first 
two moments (mean and variance [45,50]) and the autocor- 
relation function are time-invariant, albeit a stationary pro- 
cesses may be auto-correlated [60]. In the industrial/ 
engineering sphere, practitioner response to process auto- 
correlation [61] was to embrace a time-series paradigm and 
apply SPC methods to the residuals of a formal time-series 
model [45,62], albeit there were different tactical ap- 
proaches [15,63]; or to develop modified control limit 
schemes [64-66]. It is instructive to note that the non- 
model based EWMAST chart (EWMA chart for stationary 
processes [66]), recommended by Winkel and Zhang [11], 
pre-supposes a stationary (not "in-control") process. In a 
systematic review of the application of statistical process 
control in healthcare, Thor et al. [5] adduced only one lit- 
erature reference [67] and a calendar year 2003 monograph 
which discussed autocorrelation in medical series. As ar- 
gued by Alwan and Roberts [45], systematic non-random 
patterns in series make separation of the classic common 
and special causes difficult, as departures from control, 
nominally traceable to special causes, are confounded by 



autocorrelation and, in the current series, seasonality. Two 
further concerns were raised by the authors; first, the un- 
due emphasis placed upon normality and the (erroneous) 
assumption that "approximate normality" implies a state 
of statistical control; and second, in the presence of a well- 
fitting time series model with residuals consistent with 
white-noise ("randomness"), it is "...futile to search for 
departures from statistical control and their corresponding 
special causes...". The latter caution resonates with the 
current finding of frequent signalling of the raw mortality 
series compared with in-control residuals from an appo- 
site time series model; with respect to the error process, 
such signalling represents false positivity [13]. 

Cook and co-workers, "...explicitly compare [d] EWMA 
(observed) and EWMA (predicted) ...[with] thresholds 
around the EM WA (predicted)...", employing the EWMA 
(A = 0.005-0.020) to "...effectively attenuate noise in the 
data and smooth an erratic but unbiased risk model" [25], 
although no criteria of "erratic" were provided. Smoothed 
control limits for the expected series were also utilised 
in a review paper by Cook et al. ([68]) and Pilcher et al. 



EWMA /. = 0.16, putative 1SD mean increase in GARCH residuals 
Raw mortality 




Calendar time: monthly 
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Raw mortality 




^ x^, v> 



Calendar time: monthly 



EWMA /. = 0.42, putative 2SD mean increase in GARCH residuals 
Raw mortality 



-.3- 




v^, v^, v^, v^, v> V J V J V J 



Calendar time: monthly 



Figure 9 Residual-EWMA control charts (3 SE control limits) for projected 1 (upper left panel), 2 (upper right pane) and 3 (lower left 
panel) SD increase of residual mean of the GARCH model (for raw mortality series); model residuals (navy line), upper control limit 
(red line), lower control limit (green line). 
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([69], X = 0.005), albeit the data structure differed; 
sequential plotting of each patient admission versus 
monthly mortality rates in the current paper. Our focus 
and methodology were different, in that we were 
concerned to both understand and formally model the 
"noise in the data". This being said, in the current series, 
the smoothed EWMA (X = 0.51) raw series (Figure 4) was 
demonstrated to signal using 3 standard error expected 
mortality control limits. 

The sophistication of time-series modelling in stand- 
ard statistical software packages makes the formal ana- 
lyses of the current study feasible; in particular, 
automated routines for application of time series models 
[70]. However, for the application of appropriate SPC to 
mortality series from multiple ICUs in a data-base, there 
are unresolved statistical issues [71,72]. From the per- 
spectives of this study, a multivariate approach may be 
established using more conventional estimators (multi- 
variate GARCH [73] and vector autoregression models 
[74]) or by newly described hierarchical/functional time 
series [75,76]. 

Conclusions 

The underlying data generating process of monthly mor- 
tality series at the ICU level displayed autocorrelation 
and seasonality, with volatility evident in the raw mortal- 
ity series. Failure to accommodate these characteristics 
by SPC measures resulted in false-positive signalling. A 
time series approach to SPC, using residual control charts, 
would appear to resolve such issues. 

Competing interests 

The authors declare that they have no competing interests. 
Authors' contributions 

The study was conceived, designed, (data)-analysed, written and critically 
revised jointly by both authors (JLM, PJS). Both authors read and approved 
the final manuscript. 

Acknowledgements 

ANZICS Centre for Outcome and Resource Evaluation (CORE) of the Australian 
and New Zealand Intensive Care Society (ANZICS): 

Australian and New Zealand Intensive Care Society, Carlton, Victoria 3053, Australia. 
Author details 

department of Intensive Care Medicine, The Queen Elizabeth Hospital, 28 
Woodville Road, Woodville, SA 501 1, Australia. 2 School of Mathematical 
Sciences, University of Adelaide, Adelaide, SA 5000, Australia. 

Received: 7 January 2013 Accepted: 17 May 2013 
Published: 24 May 2013 

References 

1 . Montgomery DC: Quality Improvement in the Modern Business Environment. 

In Introduction ot Statistical Quality Control. 7th edition. Edited by Montgomery 
DC. Hoboken, NJ: John Wiley & Sons, Inc; 2013:3-47. 

2. Woodall WH: The Use of control charts in health-care and public-health 
surveillance. J Qual Technol 2006, 38:89-104. 

3. Benneyan JC, Lloyd RC, Plsek PE: Statistical process control as a tool for 
research and healthcare improvement. Qual Saf Health Care 2003, 12:458-464. 



4. Mohammed MA, Worthington P, Woodall WH: Plotting basic control 
charts: tutorial notes for healthcare practitioners. Qual Saf Health Care 
2008, 17:137-145. 

5. Thor J, Lundberg J, Ask J, Olsson J, Carli C, Haerenstam KP, Brommels M: 
Application of statistical process control in healthcare improvement: 
systematic review. Qual Saf Health Care 2007, 16:387-399. 

6. Collins G, Jibawi A, McCulloch P: Control chart methods for monitoring 
surgical performance: a case study from gastro-oesophageal surgery. 
Ejso 2011,37:473-480. 

7. Cook DA, Steiner SH, Cook RJ, Farewell VT, Morton AP: Monitoring the 
evolutionary process of quality: risk-adjusted charting to track outcomes 
in intensive care. Crit Care Med 2003, 31:1676-1682. 

8. Duclos A, Voirin N, Touzet S, Soardo P, Schott AM, Colin C, Peix JL, Lifante JC: 
Crude versus case-mix-adjusted control charts for safety monitoring in 
thyroid surgery. Qual Saf Health Care 2010, 19:1-4. 

9. Kirkham JJ, Bouamra O: The use of statistical process control for 
monitoring institutional performance in trauma care. J Trauma 2008, 
65:1494-1501. 

10. Rogers CA, Reeves BC, Caputo M, Ganesh JS, Bonser RS, Angelini GD: 
Control chart methods for monitoring cardiac surgical performance and 
their interpretation. J Thorac Cardiovasc Surg 2004, 128:81 1-819. 

1 1 . Winkel P, Zhang NF: Statistical development of Quality in Medicine. Chichester, 
West Sussex: John Wiley & Sons Ltd; 2007. 

12. Montgomery DC: Other univariate statistical process monitoring and 
control techniques. In Introduction ot Statistical Quality Control. 7th edition. 
Edited by Montgomery DC. Hoboken, NJ: Wiley; 2013:448-508. 

13. Alwan LC: Effects of autocorrelation on control chart performance. 
Commun Stat Theory Methods 1992, 21:1025-1049. 

14. Berthouex PM, Hunter WG, Pallesen L: Monitoring sewage-treatment 
plants - some quality-control aspects. J Qual Technol 1978, 10:139-149. 

15. Montgomery DC, Mastrangelo CM: Some statistical process-control 
methods for autocorrelated data. J Qual Technol 1991, 23:179-193. 

16. Wardell DG, Moskowitz H, Plante RD: Control charts in the presence of 
data correlation. Manag Sci 1992, 38:1084-1 105. 

17. Winkel P, Zhang NF: Serial correlation of quality control data - on the use 
of proper control charts. Scand J Clin Lab Invest 2004, 64:195-203. 

18. Winkel P, Zhang NF: Control charts for autocorrelated data. In Statistical 
development of Quality in Medicine. Edited by Winkel P, Zhang NF. 
Chichester, West Sussex: Wiley; 2007:92-1 10. 

19. Moran JL, Solomon PJ, Adult Database Management Committee (ADMC) of 
the Australian and New Zealand Intensive Care Society (ANZICS): Conventional 
and advanced time series estimation: application to the Australian and New 
Zealand Intensive Care Society (ANZICS) adult patient database, 1993-2006. 
J Eval Clin Pract 2011,17:45-60. 

20. Moran J, Solomon P: Mortality and Intensive Care volume in ventilated 
patients, 1995-2009, in the Australian and New Zealand bi-national adult 
patient intensive care database. Crit Care Med 2012, 40:800-812. 

21. Stow PJ, Hart GK, Higlett T, George C, Herkes R, McWilliam D, Bellomo R: 
Development and implementation of a high-quality clinical database: 
the Australian and New Zealand Intensive Care Society adult patient 
database. J Crit Care 2006, 21:133-141. 

22. ANZICS Centre for Outcome and Resource Evaluation (CORE) of the 
Australian and New Zealand Intensive Care Society (ANZICS): APD Data 
Dictionary: Version 3.2.1 Updated February 2012. [http://www.anzics.com.au/ 
core/data-collection-tools], Accessed June 30th 2012. 

23. Knaus WA, Wagner DP, Draper EA, Zimmerman JE, Bergner M, Bastos PG, 
Sirio CA, Murphy DJ, Lotring T, Damiano A: The APACHE III prognostic 
system. Risk prediction of hospital mortality for critically ill hospitalized 
adults. Chest 1991, 100:1619-1636. 

24. Rabe-Hesketh S, Skrondal A: Multilevel and Longitudinal Modeling Using 
Stata. 2nd edition. College Station, TX: Stata Press; 2008. 

25. Cook DA, Coory M, Webster RA: Exponentially weighted moving average 
charts to compare observed and expected values for monitoring risk- 
adjusted hospital indicators. BMJ Qual Saf 201 1, 20:469-474. 

26. Shiskin J: Decomposition of economic time series. Science 1958, 
128:1539-1546. 

27. R Development Core Team: R : A Language and Environment for Statistical 
Computing. Vienna, Austria: R Foundation for Statistical Computing; 2012 
[http://www.R-project.org]. 

28. Stoffer D: Applied Statistical Time Series Analysis ("astsa"): R package (V 1.1). 
[http://www.stat.pitt.edu/stoffer/tsa3/], Accessed 12th August 2012. 



Moran and Solomon BMC Medical Research Methodology 2013, 13:66 
http://www.biomedcentral.com/1471-2288/13/66 



Page 12 of 12 



29. Montgomery DC: Cumulative sum and exponentially weighted moving 
averaged control charts. In Introduction ot Statistical Quality Control. 7th 
edition. Edited by Montgomery DC. Hoboken, NJ: Wiley; 2013:413-447. 

30. StataCorp: lime Series Manual: Release 12. College Station, TX: StataCorp LP; 201 1 . 

31. Reynolds MR, Stoumbos ZG: Comparisons of some exponentially 
weighted moving average control charts for monitoring the process 
mean and variance. Technometrics 2006, 48:550-567. 

32. Montgomery DC: Methods and Philosophy of Statistical Process Control. 
In Introduction ot Statistical Quality Control. 7th edition. Edited by 
Montgomery DC. Hoboken, NJ: Wiley; 2013:187-233. 

33. Zhang NF: Statistical control for autocorrelated data. Proc Soc Photo Opt 
Instrum Eng 1999, 3742:65-70. 

34. StatPoint Technologies Inc: STATGRAPHICS Centurion XVI. 1. 201 2. Warrenton, 
Virginia USA. 

35. Fang Y, Zhang J: Performance of control charts for autoregressive 
conditional heteroscedastic processes. J Appl Stat 1999, 26:701-714. 

36. Schipper S, Schmid W: Control charts for GARCH processes. Nonlinear 
Anal-Theor 2001, 47:2049-2060. 

37. Sh urn way RH, Stoffer DS: ARIMA models. In Time Series Analaysis and Its 
Applications With R Examples, third edition. Edited by Shumway RH, Stoffer 
DS. New York, NY: Springer Science+Business Media, LLC; 2011:83-172. 

38. Wang S-H, Hafner C: Estimating autocorrelations in the presence of 
deterministic trends. J Time Ser Econom 201 1, 3. Articled http://www. 
degruyter.com/view/j7jtse.201 1 .3.2/jtse.201 1 .3.2.1 022/jtse.201 1 .3.2.1 022.xml. 

39. Nelson CR, Plosser CR: Trends and random walks in macroeconmic time 
series : some evidence and implications. J Monet Econ 1982, 10:139-162. 

40. Pierce DA: Trend and autocorrelation. Commun Stat 1975, 4:163-1 75. 

41. Shumway RH, Stoffer DS: Additional time domain topics. In Time Series 
Analaysis and Its Applications With R Examples, third edition. Edited by 
Shumway RH, Stoffer DS. New York, NY: Springer Science+Business Media, 
LLC; 2011:267-318. 

42. Tsay RS: Asset volatility and Volatility models. In An Introduction to Analysis of 
Financial Data with R. Edited by Tsay RS. Hoboken, NJ: Wiley; 2013:176-241. 

43. Karlsson S: ARMADIAG: Stata module to compute post-estimation residual 
diagnostics for time series; 2009. [http://econpapers.repec.org/scripts/search. 
asp?ft=armadiag], Accessed May 2009. 

44. Engle R: GARCH 101: the use of ARCH/GARCH models in applied 
econometrics. J Econ Perspect 2001, 15:157-168. 

45. Alwan LC, Roberts HV: Time-series modeling for statistical process-control. 
J Bus Econ Stat 1988, 6:87-95. 

46. Koehler AB, Marks NB, O'Connell RT: EWMA control charts for autoregressive 
processes. J Oper Res Soc 2001, 52:699-707. 

47. Tsay RS: An Introduction to Analysis of Financial Data with R. John Wiley & 
Sons, Inc: Hoboken, NJ; 2013. 

48. Kuha J: AIC and BIC: comparisons of assumptions and performance. 
Sociol Method Res 2004, 33:188-229. 

49. StataCorp: tssmooth exponential -Single-exponential smoothing. In Time 
Series Manual: Release 12. College Station, TX; 201 1:477-484. 

50. Alwan LC, Roberts HV: The problem of misplaced control limits. Appl Stat-J 
/toy St C 1995, 44:269-278. 

51. Bisgaard S, Kulahci M: Quality quandaries: the effect of autocorrelation on 
statistical process procedures. Qual Eng 2005, 17:481-489. 

52. Box G, Narasimhan S: Rethinking statistics for quality control. Qual Eng 
2010, 22:60-72. 

53. Apley DW, Lee HQ Robustness comparison of exponentially weighted 
moving-average charts on autocorrelated data and on residuals. J Qual 
Technol 2008, 40:428-447. 

54. Lu CW, Reynolds MR: EWMA control charts for monitoring the mean of 
autocorrelated processes. J Qual Technol 1999, 31:166-188. 

55. Reynolds MR, Lu CW: Control charts for monitoring processes with 
autocorrelated data. Nonlinear Anal-Theor 1997, 30:4059-4067. 

56. Human SW, Kritzinger P, Chakraborti S: Robustness of the EWMA control 
chart for individual observations. J Appl Stat 201 1, 38:2071-2087. 

57. Vasilopoulis AV, Stamboulis AP: Modification of control chart limits in the 
presence of data correlation. J Qual Technol 1 978, 1 0:20-30. 

58. Shahian DM, Williamson WA, Svensson LG, Restuccia JD, DAgostino RS: 
Applications of statistical quality control to cardiac surgery. Ann Thorac 
Surg 1996, 62:1351-1358. 

59. Tennant R, Mohammed MA, Coleman JJ, Martin U: Monitoring patients 
using control charts: a systematic review. Int J Qual Health Care 2007, 
19:187-194. 



60. Kirchgassner G, Wolters J: Introduction to Modern Times Series Analysis. Berlin: 
Springer; 2008. 

61. Lwin T: Parameter estimation in first-order autoregressive model for 
statistical process monitoring in the presence of data autocorrelation. 

J Stat Plan Infer 2011, 141:2556-2575. 

62. Roberts HV, Tsay RS: Making control charts more effective by time series 
analysis: three illustrative applications. Commun Stat-Theory 1996, 
25:2767-2796. 

63. Runger GC: Assignable causes and auto correlation: control charts for 
observations or residuals? J Qual Technol 2002, 34:165-170. 

64. Apley DW: Time series control charts in the presence of model uncertainty. 
J ManufSci E-TAsme 2002, 1 24:891 -898. 

65. Lee HC, Apley DW: Improved design of robust exponentially weighted 
moving average control charts for autocorrelated processes. Qual Reliab 
Engng Int 2011, 27:337-352. 

66. Zhang NF: A statistical control chart for stationary process data. 
Technometrics 1 998, 40:24-38. 

67. Solodky C, Chen HG, Jones PK, Katcher W, Neuhauser D: Patients as partners in 
clinical research - a proposal for applying quality improvement methods to 
patient care. Med Care 1 998, 36:AS1 3-AS20. 

68. Cook DA, Duke G, Hart GK, Pilcher D, Mullany D: Review of the application 
of risk-adjusted charts to analyse mortality outcomes in critical care. 
Crit Care Resusc 2008, 10:239-251. 

69. Pilcher DV, Hoffman T, Thomas C, Ernest D, Hart GK: Risk-adjusted 
continuous outcome monitoring with an EWMA chart: could it have 
detected excess mortality among intensive care patients at Bundaberg 
Base Hospital? Crit Care Resusc 2010, 12:36-41. 

70. Hyndman RJ, Khandakar Y: Automatic time series forecasting: the forecast 
package for R. J Stat Softw 2008, 27:1 -22. 

71 . Bottle A, Aylin P: Predicting the false alarm rate in multi-institution mortality 
monitoring. J Oper Res Soc 201 1 , 62:1 71 1 -1 71 8. 

72. Marshall T, Mohammed MA, Rouse A: A randomized controlled trial of 
league tables and control charts as aids to health service decision-making. 
Int J Qual Health Care 2004, 16:309-315. 

73. Bauwens L, Laurent S, Rombouts JV: Multivariate GARCH models: a survey. 
J Appl Econ 2006,21:79-109. 

74. Pan X, Jarrett JE: Why and how to use vector autoregressive models for 
quality control: the guideline and procedures. Qual Quant 2012, 46:935-948. 

75. de Silva A, Hyndman RJ, Snyder R: The vector innovations structural time 
series framework: a simple approach to multivariate forecasting. Stat 
Model 2010, 10:353-374. 

76. Hyndman RJ, Ahmed RA, Athanasopoulos G, Shang HL: Optimal combination 
forecasts for hierarchical time series. ComputStat Data An 201 1, 55:2579-2589. 



doi:1 0.1 1 86/1 471 -2288-1 3-66 

Cite this article as: Moran and Solomon: Statistical process control of 
mortality series in the Australian and New Zealand Intensive Care Society 
(ANZICS) adult patient database: implications of the data generating 
process. BMC Medical Research Methodology 201 3 1 3:66. 



Submit your next manuscript to BioMed Central 
and take full advantage of: 

• Convenient online submission 

• Thorough peer review 

• No space constraints or color figure charges 

• Immediate publication on acceptance 

• Inclusion in PubMed, CAS, Scopus and Google Scholar 

• Research which is freely available for redistribution 



Submit your manuscript at 
www.biomedcentral.com/submit 



o 



BioMed Central 



